function [outq,outm,outcell]=dsge_spectrum_mixfreq(GG,RR,ZZ,SDX,specrep,stanames,shonames,...
    obsdesc,lowper,highper,npoints,delta,outpath,figname) 

% function [outq,outcell,hand_vec]=dsge_spectrum_mixfreq(GG,RR,ZZ,SDX,specrep,stanames,shonames,...
% obsdesc,lowper,highper,npoints,delta,outpath,figname) 
%
%% DSGE_SPECTRUM_MIXFREQ
% 
%% 1. Purpose
% Compute the spectrum of a DSGE whose solution is in (time
% unit,frequency,periodicity) given by $$ (t,\omega,p_{m}) $$ following: 
% 
% A. Smoothing through the filter $$ B(L) $$ 
%
% B. Changing to the triplet $$ (\tau,\lambda,p_{h}) $$ 
%    where $$ \tau $$ is the lower, e.g. sampling time unit and where there
%    are $$ \delta $$ t-periods per every one $$ \tau $$ period. 
%    Put differently, the sampling unit is for $$ \tau=1,...,T $$ with $$
%    \tau=\delta t $$ and $$ t=1,....,T \delta $$. 
%
% In principle, the filter $$ B(L) $$ can be either an average or a sum of lags,
%  depending on the value of the variable *CASE_TAG* for the observables
%  and *CASE_TAGST* for the states. If equal to 2, then the filter is an
%  average of $$ \delta $$ periods while if equal to 1, just a sum. 
% 
%  * NOTE: Currently assume AVERAGE, i.e. CASE_TAG =2  for all series *
% 
%% 2. Inputs 
% 
% GG,RR,ZZ amd SDX are the usual solution matrices of the DSGE 
% 
% Note: this will depend on whether the model is GTI or GTV. 
% 
% i)For the GTV, ZZ this must include a mapping to the states that are part of
% the cumulator, e.g. ZZ=[Z;A] 
% 
% ii) For the GTI, ZZ=Z. 
% 
% *SPECREP*: cell, with names of states wish to report. All observables are reported
% automatically. If ' lev' then it will be multiplied by the (1-L)^(-1) filter. 
% 
% *STANAMES*: cell with the names of all states. 
%
% *SHONAMES*: cell with the names of all shocks. 
% 
% *OBSDESC*: cell with the names of observables, in the ZZ matrix. 
% 
% *LOWPER*: scalar, lower periodicity > 2, for the QUARTERLY model. 
%
% *HIGHPER*: scalar, higher periodicity, for the QUARTERLY model.
%
% *NPOINTS*: scalar, number of points used in the approximation, 1000 should provide
%           a reasonable approximation. 
% 
% *DELTA*: number of periods in the short, i.e. model, frequency, at which
% some data are observed, i.e. sampling frequency. Should be 3 for M to Q. 
%
% *OUTPATH*: if provided, a table will be created and stored in . 
% 
%% 3. Output in Structure OUTQ
% 
% *3.I Observables* 
%
% *SPEC_QOBS*= [NZ NZ NPOINTS] gives $$ S_{y}(\lambda_{j}) $$ 
%
% *SPECDEC_QOBS* =[NZ NX NPOINTS]  gives $$ S_{y_i,y_i}(\lambda_{j}) $$ 
% i.e, every row is a decomposition across all shocks of the variance at that frequency.  
% 
% *SPECVFBDEC_QOBS* =[NZ 1] gives $$ \sum S_{y}(\lambda_{j})*Area_{j} $$
% 
% *SPECVFBDEC_QOBS* =[NZ NX] each row gives the ratio 
% $$ \sum S^j_{y_i,y_i}(\lambda_{j})*Area_{j} \div 
%     \sum S_{y_i,y_i}(\lambda_{j})*Area_{j } $$ 
% 
% *3.II States Unstranformed*
%
% Following are the same as their _QOBS counterparts, but for the states
% one wishes to report
%
% *SPEC_QST*= [NSTQ NSTQ  NPOINTS], see SPEC_QOBS
%
% *SPECDEC_QST* =[NSTQ NX NPOINTS], see SPECDEC_QOB
%
% *SPECVFBDEC_QST* =[NSTQ 1], see SPECVFBDEC_QOBS
% 
% *SPECVFBDEC_QST* =[NSTQ  NX] see SPECVFBDEC_QOBS
%
% *3.III States in Levels* 
% 
% *SPEC_QSTLEV*= [NSTLEV NSTLEV NPOINTS], see SPEC_QOBS
%
% *SPECDEC_QSTLEV* =[NSTLEV NX NPOINTS], see SPECDEC_QOB
%
% *SPECVFBDEC_QSTLEV* =[NSTLEV 1], see SPECVFBDEC_QOBS
% 
% *SPECVFBDEC_QSTLEV* =[NSTLEV NX] see SPECVFBDEC_QOBS
%  
%% 4. Notes and Related Codes 
% 
% Related codes: VARDECOM_SPEC.m and DSGE_SPECTRUM for models without TAG 
%
% Created April 9 2010 
% Last Modified: 
% Alejandro Justiniano April 28 2010 
%

%% I. Indicators 
% 
% Obtained Indicators of which rows of the state vector to report and
% whether to transform to a level or not. 
%
% Dimensions:
% *NSTQ*:   Number of states to report 
% *NSTLEV*: Number of states to report in levels 
% *POSLEVIND*: Indicator of the position WITHIN the NSTQ states, to
% aggregate to a level 
%
% STSPEC        Position in STANAMES of the states to report 
%
% SPECLEVIND    [0 1] Indicator of which of these states will be in levels
%
% STNAMERAW     Is the name of the RAW, i.e. untransformed states 
%
% 
cucd=pwd; 
if nargin < 14 || isempty( figname ) 
    flag_figure=0; 
else 
    flag_figure=1;
end 
    
if nargin < 13 || isempty(outpath); 
    flag_save=0; 
else 
    flag_save=1; 
end 

stanames=stanames(:); 
shonames=shonames(:);

if ~isempty(specrep)
    flag_repstates=1;
    nstq     =length(specrep);
    stspec    =zeros(nstq,1);
    speclevind=zeros(nstq,1);
    stnameraw =cell(nstq,1);
    for ii=1:nstq;
        tempos=findstr(' lev',specrep{ii} );
        if ~isempty( tempos);
            tempstr=specrep{ii};
            stnameraw(ii)={tempstr(1:tempos-1)};
            stspec(ii) =cellposition( stnameraw(ii) , stanames );
            speclevind(ii)=1;
        else
            stspec(ii) =cellposition( specrep(ii), stanames );
            speclevind(ii)=0;
            stnameraw(ii)=specrep(ii);
        end
    end
    clear tempos tempstr;
    poslevind=find( speclevind == 1);
    nstlev=length(poslevind);
    
    disp('TAG Spectrum: States to report and level indicator')
    printcell( [stanames(stspec) num2cprec(speclevind)] )
    
else
    flag_repstates=0;
    disp('No states will be reported')
    stspec=[];
    
end
clear temp;
    
%% II. Quarterly and Monthly periodicities and frequencies 
% 
% Given the inputs $$ (p_{qL},p_{qH}) $$ corresponding to the quarterly Low
% and High periodicities, find: 
%
% a. A grid of periodicities $$ p_{q,j}  j=1,...,NPOINTS $$ 
% 
% b. The associated quarterly frequencies $$ \lambda_{j}=2 \pi / p_{q,j} $$ 
%
% c. The associated monthly   frequencies $$  \omega_{j,h}= ( \lambda_{j} +
% 2h \pi )/\delta $$ for $$ h=0,...,\delta-1 $$
% 
% Note: for diagonal elements only, this formula is identical to the one
% provided by Granger and Siklos, JE 1995, for whom the frequencies that
% define the time aggregated spectrum correspond to 
% $$  \omega_{j,0}= ( \lambda_{j} / delta ) ; \omega_{j,1}= ( \lambda_{j} +
% 2 \pi )/\delta ; \omega_{j,2}= ( \lambda_{j} -
% 2 \pi )/\delta  $$
% 
% This is because the diagonal elements of the multivariate spectrum at $$ \omega_{j,2}= ( \lambda_{j} -
% 2 \pi )/\delta $$ equal those at $$  ( \lambda_{j} -
% 2 \pi )/\delta + 2 \pi = ( \lambda_{j} +
% 2 (\delta-1) \pi )/\delta $$ as in Harvey's original formula. 
% 
% Therefore, the computations provided here *are ONLY valid for the diagonal elements of the spectrum* 
% 
% 

%% III. Grid for $$ (\lambda , \omega) $$ 
% 
% Note: In the case w/o time aggregation, take N+2 points and then truncate to the interior N points.  
% In this way, the AREAS of the basis of the rectangles are given by 
% [W(1)-W(0), W(2)-W(1),.....,W(N)-W(N-1)] 
% where W(0)  =2pi/HIGHPER 
%       W(N+1)=2pi/LOWPER , i.e. W(N) is strictly an interior point and the
%       area between W(N) and 2pi/LOWPER is not covered 
% We will assume the same here. This creates a bit of an ambiguity
% regarding the withd of the edges for the monthly frequency. But with 
% DELTA*NPOINTS, its influence should be very small in the monthly spectrum
% and will be ZERO for the quarterly spectrun, since the computation of the
% areas is done below, with the correct withds W(0) and W(N) defined as
% above. 
%
pqvec=linspace(lowper,highper,npoints+2);
lambda=fliplr((2*pi)./pqvec);
lambda=lambda(2:end-1);
pqvec=(2*pi)./lambda; 

% COMPUTE the OUT structure 
% -------------------------------------------------------------------------
disp('Spectral Analysis for TAG model ');
dispaj('For periodicities ',['(',num2str(lowper),',',num2str(highper),')']); 


% Follow Harvey's formula
omega=zeros(npoints,delta);
for ii=1:npoints  
    omega(ii,:)= ( lambda(ii)*ones(1,delta) + (0:1:delta-1)*2*pi )/delta; 
end 
% Obtain vector of UNIQUE frequencies for which to evaluate the spectrum 
% Note: Could be reduced even further by limiting to the [0,pi] interval,
% which would require adjusting the position indicators derived next. 
omegain=omega'; 

omegain=sort(unique( omegain(:) ),'ascend'); 
% OMEGAPOSIN gives the position of the corresponding OMEGA frequencies in
% the vector of unique frequencies OMEGAPOSIN passed for the computation of
% the spectrum 
omegaposin=zeros(npoints,delta);
for ii=1:npoints
    for jj=1:delta
        omegaposin(ii,jj)=find( omegain==omega(ii,jj) );
    end
end

%% IV. Compute Monthly Spectrum $$ S_{x,x}(\omega) $$ 
% *OUTM*: Output of the monthly spectrum 
% 
% Note: as in the case above, the computations are carried out in the $$
% \omega \in [0,\pi] $$ interval. Therefore, for the off-diagonal elements,
% should not take $$ S_{\omega}  \omega  \in [- \pi ,0] $$ to equal $$ S_{\omega + \pi} $$ since the imaginary 
% elements of the COSPECTRUM are the negative of one another. 
% All calculations are valid either for integration over a frequency band,
% for which the COSPECTRUM integrates to zero in a symmetric interval, as
% well as for the diagonal elements. 
% 
% EDGEM: obtain distance to edges. 
% Without temporal aggregation this is obtained as 
% wlb  =2*pi/highper;wub=2*pi/lowper; 
% edgew=[(freq(1)-wlb) (wub-freq(end))];
% Here, as a simplification, we will assume it is
% symmetric to median distance in the set of unique frequencies. 
% As explained above, this has no effect on the Variances of the Aggregated
% spectrum, since these are done below. 
edgewm=median( omegain(2:end)-omegain(1:end-1) )*ones(2,1); 
if any( edgewm < 0 ); error('EDGEW cannot be negative'); end 

% VARDECOM_SPEC computes Monthly Spectrum 
outm=vardecom_spec(GG,RR,SDX,ZZ,omegain,edgewm,stspec,zeros(length(stspec),1)); 

%% IV.1. Transform $$ S_{x,x}(\omega) $$ into $$ S_{y,y}(\lambda) $$
% 
% Transformation 1: Spectrum for AVERAGES of monthly series. 
%
% Compute Monthly Spectrum $$ S_{xb,xb}(\omega) $$ by multiplying through the frequency response 
% function of the filter $$ B(L) $$. It is assumed the transformation is
% the same for ALL series. 
% 
% Transformation 2: SUM the Spectrum for AVERAGES through the Folding
% Formula. I.e., change the time unit from model to sampling. 
%

%% IV.2. Observables 
%
% *SPEC_QOBS*= [NZ NZ NPOINTS] gives $$ S_{y}(\lambda_{j}) $$ 
%
% *SPECDEC_QOBS* =[NZ NX NPOINTS]  gives $$ S_{y_i,y_i}(\lambda_{j}) $$ 
% i.e, every row is a decomposition across all shocks of the variance at that frequency.  
% 
% *SPECVFBDEC_QOBS* =[NZ 1] gives $$ \sum S_{y}(\lambda_{j})*Area_{j} $$
% 
% *SPECVFBDEC_QOBS* =[NZ NX] each row gives the ratio 
% $$ \sum S^j_{y_i,y_i}(\lambda_{j})*Area_{j} \div 
%     \sum S_{y_i,y_i}(\lambda_{j})*Area_{j } $$ 
% 
nz=size(ZZ,1); 
nx=size(SDX,1); 
% STQ is the number of states that are being reported 
outq.spec_qobs      =zeros(nz,nz,length(lambda)); 
outq.specdec_qobs   =zeros(nz,nx,length(lambda)); 
outq.specvfbdec_qobs=zeros(nz,nx); 
outq.specvfb_qobs   =zeros(nz,1 ); 

%% IV.3. States 
%
% IV.3.i States Unstranformed
%
% Following are the same as their _QOBS counterparts, but for the states
% one wishes to report
%
% *SPEC_QST*= [NSTQ NSTQ  NPOINTS], see SPEC_QOBS
%
% *SPECDEC_QST* =[NSTQ NX NPOINTS], see SPECDEC_QOB
%
% *SPECVFBDEC_QST* =[NSTQ 1], see SPECVFBDEC_QOBS
% 
% *SPECVFBDEC_QST* =[NSTQ  NX] see SPECVFBDEC_QOBS
%
% IV.3.i States in Levels 
% 
% *SPEC_QSTLEV*= [NSTLEV NSTLEV NPOINTS], see SPEC_QOBS
%
% *SPECDEC_QSTLEV* =[NSTLEV NX NPOINTS], see SPECDEC_QOB
%
% *SPECVFBDEC_QSTLEV* =[NSTLEV 1], see SPECVFBDEC_QOBS
% 
% *SPECVFBDEC_QSTLEV* =[NSTLEV NX] see SPECVFBDEC_QOBS
%  
if flag_repstates==1
    outq.spec_qst      =zeros(nstq,nstq,length(lambda));
    outq.specdec_qst   =zeros(nstq,nx,length(lambda));
    outq.specvfbdec_qst=zeros(nstq,nx);
    outq.specvfb_qst   =zeros(nstq,1 );
    if nstlev > 0
        outq.spec_qstlev      =zeros(nstlev,nstlev,length(lambda));
        outq.specdec_qstlev   =zeros(nstlev,nx,length(lambda));
        outq.specvfbdec_qstlev=zeros(nstlev,nx);
        outq.specvfb_qstlev   =zeros(nstlev,1 );
    end
end
    
% EDGES 
% Recall EDGEQ=[W(1)-W(0), W(2)-W(1),.....,W(N)-W(N-1)] 
% where W(0)  =2pi/HIGHPER 
%       W(N+1)=2pi/LOWPER , i.e. W(N) is strictly an interior point and the
%       area between W(N) and 2pi/LOWPER is not covered 
edgeq=zeros(length(lambda),1); 
edgeq(1      )=lambda(1)-2*pi/highper;
edgeq(2:end  )=lambda(2:end)-lambda(1:end-1); 
%edgeq(end    )=2*pi/lowper-lambda(end); 

for ii=1:length(lambda)
        
    temp_spec_xbobs   =zeros(nz,nz,delta); 
    temp_specdec_xbobs=zeros(nz,nx,delta); 
    
    if flag_repstates==1
        temp_spec_xbst   =zeros(nstq,nstq,delta);
        temp_specdec_xbst=zeros(nstq,nx  ,delta);
    end 
    
    %%  A. Frequency Response function for an average filter
    % $$ E1= \omega \otimes 1' $$
    %
    % $$ E2= 1' \otimes [-0i -1i   -i \delta] $$
    %
    % $$ E =B(exp(-1i \omega_{h}) )= exp( E1 .* E2 ) $$
    temp_A=reshape( kron( omega(ii,:), ones(1,delta) ) ,[delta delta]) ;
    temp_B=reshape( kron( ones(1,delta), -1i*(0: (delta-1) ) ) ,[delta delta] );
    temp_C=sum( exp( temp_A.*temp_B ) , 1 )/delta;
    for jj=1:delta;
        % Position indicator of where to extract from the MONTHLY SPECTRUM
        pos =omegaposin(ii,jj);
        temp_spec_xbobs(:,:,jj)   =(abs(temp_C(jj))^2)*outm.specz(:,:,pos);
        % .SPECDECZ is the SHARE of the diagonal, hence multiply by the
        % diagonal to get the levels. VARDECOM_SPEC already check these
        % shares add up to one, but will also be checked below
        temp_specdec_xbobs(:,:,jj)=(abs(temp_C(jj))^2)*( outm.specdecz(:,:,pos).*repmat(diag(outm.specz(:,:,pos)),[1 nx]) );
        
        if flag_repstates==1
            temp_spec_xbst(:,:,jj)=(abs(temp_C(jj))^2)*outm.specst(:,:,pos);
            temp_specdec_xbst(:,:,jj)=(abs(temp_C(jj))^2)*( outm.specdecst(:,:,pos).*repmat(diag(outm.specst(:,:,pos)),[1 nx]) );
        end
                
    end
    %% B. FOLDING FORMULA 
    % $$ S_y(\lambda) = \sum{ |B( exp( -i \omega_{h} ) )|^2
    % S_x( \omega_{h} )   } \div \delta $$ 
    outq.spec_qobs(:,:,ii)       =reshape( sum(temp_spec_xbobs,3   ) , [nz nz] )/delta; 
    outq.specdec_qobs(:,:,ii)    =reshape( sum(temp_specdec_xbobs,3) , [nz nx] )/delta; 

    if flag_repstates==1
        outq.spec_qst(:,:,ii)     =reshape( sum(temp_spec_xbst,3   ) , [nstq nstq] )/delta;
        outq.specdec_qst(:,:,ii)  =reshape( sum(temp_specdec_xbst,3) , [nstq nx  ] )/delta;
        
        if nstlev > 1 
            %% C. Frequency reponse for $$ [ 1-exp(-i \lambda_{j}  ]^(-1)
            % frfinv=(1/( 2-2*cos(lambda(ii)) ));             
            outq.spec_qstlev(:,:,ii)   =(1/( 2-2*cos(lambda(ii)) ))*outq.spec_qst(poslevind,poslevind,ii);
            outq.specdec_qstlev(:,:,ii)=(1/( 2-2*cos(lambda(ii)) ))*outq.specdec_qst(poslevind,:,     ii);
        end 
        
    end
    
    % =====================================
    % Check decompositions are correct
    % ======================================
    maxdif=max( abs( squeeze( sum( outq.specdec_qobs(:,:,ii) , 2 ) ) - diag( outq.spec_qobs(:,:,ii) ) ) );
    if maxdif > 1e-5
        error('Could not decompose the Spectrum across frequencies for the quarterly aggregated OBSERVABLES')
    end    
    if flag_repstates==1
        maxdif=max( abs( squeeze( sum( outq.specdec_qst(:,:,ii) , 2 ) ) - diag( outq.spec_qst(:,:,ii) ) ) );
        if maxdif > 1e-5
            error('Could not decompose the Spectrum across frequencies for the quarterly aggregated STATES')
        end        
        if nstlev > 1            
            maxdif=max( abs( squeeze( sum( outq.specdec_qstlev(:,:,ii) , 2 ) ) - diag( outq.spec_qstlev(:,:,ii) ) ) );
            if maxdif > 1e-5
                error('Could not decompose the Spectrum across frequencies for the quarterly aggregated LEVEL STATES ')
            end
            
        end
    end            
    
    %% D. VARIANCE FREQUENCY BAND (VFB) and it's Decomposition 
    outq.specvfb_qobs   =outq.specvfb_qobs+2*diag(squeeze(outq.spec_qobs(:,:,ii)))*edgeq(ii); 
    outq.specvfbdec_qobs=outq.specvfbdec_qobs+2*outq.specdec_qobs(:,:,ii)*edgeq(ii); 

    if flag_repstates==1
        
        outq.specvfb_qst  =outq.specvfb_qst+2*diag(squeeze(outq.spec_qst(:,:,ii)))*edgeq(ii);
        outq.specvfbdec_qst=outq.specvfbdec_qst+2*outq.specdec_qst(:,:,ii)*edgeq(ii);
        
        if nstlev > 1
            
            % frfinv=(1/( 2-2*cos(lambda(ii)) ));
            outq.specvfb_qstlev  =outq.specvfb_qstlev+2*diag(squeeze(outq.spec_qstlev(:,:,ii)))*edgeq(ii);
            outq.specvfbdec_qstlev=outq.specvfbdec_qstlev+2*outq.specdec_qstlev(:,:,ii)*edgeq(ii);
        end
        
    end

    
    
end 
% =====================================================================
% Check total variance decomposition is correct 
maxdif=max( abs( sum( outq.specvfbdec_qobs, 2 )  - outq.specvfb_qobs ) );
if maxdif > 1e-5
    error('Could not decompose the Variance Frequency Band for the quarterly aggregated observables')
end
outq.specvfbdec_qobs=outq.specvfbdec_qobs./repmat( outq.specvfb_qobs , [1 nx] ); 
if max( abs( sum(outq.specvfbdec_qobs,2)-1) > 1e-5 ) 
    error('Could not decompose the Variance Frequency Band for the quarterly aggregated observables')
end 

if flag_repstates==1
    maxdif=max( abs( sum( outq.specvfbdec_qst, 2 )  - outq.specvfb_qst ) );
    if maxdif > 1e-5
        error('Could not decompose the Variance Frequency Band for the quarterly aggregated STATES ')
    end
    outq.specvfbdec_qst=outq.specvfbdec_qst./repmat( outq.specvfb_qst , [1 nx] );
    if max( abs( sum(outq.specvfbdec_qst,2)-1) > 1e-5 )
        error('Could not decompose the Variance Frequency Band for the quarterly aggregated STATES')
    end
    if nstlev > 1 
        maxdif=max( abs( sum( outq.specvfbdec_qstlev, 2 )  - outq.specvfb_qstlev) );
        if maxdif > 1e-5
            error('Could not decompose the Variance Frequency Band for the quarterly aggregated STATES ')
        end
        outq.specvfbdec_qstlev=outq.specvfbdec_qstlev./repmat( outq.specvfb_qstlev , [1 nx] );
        if max( abs( sum(outq.specvfbdec_qstlev,2)-1) > 1e-5 )
            error('Could not decompose the Variance Frequency Band for the quarterly aggregated STATES')
        end        
    end     
end 



disp('Finished Spectral Analysis' ); 

%% V. Output table 
%  
% OUTCELL shows the average decomposition over the frequency band 
% 
%empstr=['(',num2str(lowper),',',num2str(highper),')']; 
outcell=emptycell(4,6); 
outcell(1,1         )={extrsubf(outpath,cd)}; 
outcell(2,1:4       )={'Aggregated','Sums','DELTA (INT)',num2str(delta)}; 
outcell(3,[1:2  4:5])={'Lower Per',num2str(lowper),'npoints',num2str(npoints)}; 
outcell(4, 1:2      )={'High  Per',num2str(highper)                          };
try
    outspec=crtcell( outq.specvfbdec_qstlev ,  specrep(poslevind) , shonames ,{'Spec Dec LEVELS'} );
    outcell=merge_cells( outcell, outspec , 1 );
catch
    disp(['No Levels involved'])
end
outspec=crtcell( outq.specvfbdec_qst   ,  stnameraw          , shonames ,{'Spec Dec ST Select'} ); 
outcell=merge_cells( outcell, outspec , 1 ); 
outspec=crtcell( outq.specvfbdec_qobs      ,  obsdesc            , shonames ,{'Spec Dec Z Raw'} ); 
outcell=merge_cells( outcell, outspec , 1 ); 

if flag_save==1
    cd(outpath)
    xlswrite('variance decomposition tab',outcell,'spectrum');
    cd(cucd); 
end
